########################################################################
########################################################################
# Replication code for:
#
# @article{Zhukov2016JCE,
# title={Trading hard hats for combat helmets: The economics of rebellion in eastern Ukraine},
# author={Zhukov, Yuri M},
# journal={Journal of Comparative Economics},
# volume={44},
# number={1},
# pages={1--15},
# year={2016},
# publisher={Elsevier}
# }
#
# Comments, questions: zhukov-at-umich-dot-edu
########################################################################
########################################################################

rm(list=ls())
setwd("~/Dropbox/REPBIAS/JCE_Replication/")
source("functions.R")

# Load maps
load("Data/UKR_adm2_Donbass.RData")
load("Data/UA_PPL_v2.RData")
gn <- gn[order(as.numeric(as.character(gn$geonameid))),]


#######################
# Summary stats
#######################

#####
# Municipality-level
#####

# Correlation matrix
lv <- "_v"
load(paste0("Data/Combined_1PD",lv,".RData"))
data <- data.v; rm(data.v)
ind <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
eth1 <- c("RUS_CENSUS2001")
ter <- c("FOREST","dem","DIST2ROAD")
dem <- c("POP")
rus <- c("DIST2RUS")
pol <- c("V2010A_YANUKOVYCH")
X.varz <- c(ind,eth1,ter,dem,rus,pol)
X.labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010")
#png(file=paste0("Output/cormat",lv,".png"),width=4,height=4,units="in",res=300)
cormat(data,vars=X.varz,labels=X.labz)
#dev.off()

# Summary stats table
load("Data/SurvMat_v.RData")
data <- merge(data,surv.mat,by="geonameid")
y.varz <- c("REB_ALL_b","SURV_CTR","SURV_CTR_C","SURV_LIB2","SURV_LIB2_C")
y.labz <- c("Any rebel violence","Duration until rebel control","Duration until rebel control (censoring dummy)","Duration until loss of rebel control","Duration until loss of rebel control (censoring dummy)")
X.labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010")
xtable(sumstat(data,varz=c(y.varz,X.varz),labs=c(y.labz,X.labz)),digits=2)


#####
# Municipality-week level
#####

# Correlation matrix

lv <- "_vw"
load(paste0("Data/Combined_1PD",lv,".RData"))
data <- data.vw; rm(data.vw)
ind <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
eth1 <- c("RUS_CENSUS2001")
ter <- c("FOREST","dem","DIST2ROAD")
dem <- c("POP")
rus <- c("DIST2RUS")
pol <- c("V2010A_YANUKOVYCH")
X.varz <- c(ind,eth1,ter,dem,rus,pol)
X.labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010")
#png(file=paste0("Output/cormat",lv,".png"),width=4,height=4,units="in",res=300)
cormat(data,vars=X.varz,labels=X.labz)
#dev.off()

# Summary stats table
y.varz <- c("REB_ALL")
y.labz <- c("Intensity of rebel violence")
X.labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010")
xtable(sumstat(data,varz=c(y.varz,X.varz),labs=c(y.labz,X.labz)),digits=2)


########################################################################
########################################################################
## Bayesian Model Averaging
########################################################################
########################################################################

#####################################
# Incidence of any rebel violence
# (Figure 4a)
#####################################

# Load data
lv <- "_v"
load(paste0("Data/Combined_1PD",lv,".RData"))
data <- data.v; rm(data.v)
tvar <- "WID"
data <- data[order(data$geonameid),]
disag <- data[,"geonameid"]

# Rescale variables
data$POP <- data$POP/max(data$POP,na.rm=T)
data$dem <- data$dem/max(data$dem,na.rm=T)
data$DIST2ROAD <- data$DIST2ROAD/max(data$DIST2ROAD,na.rm=T)
data$DIST2RUS <- data$DIST2RUS/max(data$DIST2RUS,na.rm=T)
data$RUS_CENSUS2001 <- data$RUS_CENSUS2001/100
data$V2010A_YANUKOVYCH <- (data[,"V2010A_YANUKOVYCH"]-min(data[,"V2010A_YANUKOVYCH"]))/(max(data[,"V2010A_YANUKOVYCH"])-min(data[,"V2010A_YANUKOVYCH"]))

# Create interactions
data$RUS_MAJ <- 1*(data$RUS_CENSUS2001>.5)
data$MCH_RUS <- data$IndustryMachinery_Emp_PROP*data$RUS_CENSUS2001
data$MIN_RUS <- data$IndustryMining_Emp_PROP*data$RUS_CENSUS2001
data$MTL_RUS <- data$IndustryMetals_Emp_PROP*data$RUS_CENSUS2001 
data$MCH_RUS_b <- data$IndustryMachinery_Emp_PROP*data$RUS_MAJ
data$MIN_RUS_b <- data$IndustryMining_Emp_PROP*data$RUS_MAJ
data$MTL_RUS_b <- data$IndustryMetals_Emp_PROP*data$RUS_MAJ 

# Choose covariates
ind <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
eth1 <- c("RUS_CENSUS2001")
eth1_b <- c("RUS_MAJ")
int <- c("MCH_RUS","MIN_RUS","MTL_RUS")
int_b <- c("MCH_RUS_b","MIN_RUS_b","MTL_RUS_b")
ter <- c("FOREST","dem","DIST2ROAD")
dem <- c("POP")
rus <- c("DIST2RUS")
pol <- c("V2010A_YANUKOVYCH")

# BMA
X.varz <- c(ind,eth1_b,ter,dem,rus,pol,int_b)
y.varz <- c("REB_ALL_b")
X <- na.omit(data[,X.varz])
y <- data[row.names(data)%in%row.names(X),y.varz]
mod <- bic.glm(x=X,y=y,glm.family="binomial",strict=F,OR=10^3)
summary(mod,n.models=1,digits=1)
invlink <- mod$linkinv

# Size of model space (excluding incomplete interactions)
2^6*(2^3-3)*(2^3-4)*(2^3-4)

# Fix model weights
ix <- list(c("MCH_RUS_b","IndustryMachinery_Emp_PROP","RUS_MAJ"),c("MIN_RUS_b","IndustryMining_Emp_PROP","RUS_MAJ"),c("MTL_RUS_b","IndustryMetals_Emp_PROP","RUS_MAJ"))
mod <- bma.int(mod=mod,ix=ix)

# Plot
labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010","Machine x Russian","Mining x Russian","Metals x Russian")
#png(file=paste0("Output/BMA_mod1_inter",lv,".png"),width=4.5,height=4.5,units="in",res=300)
BMA.plot(mod,labz,dvlab="Any rebel violence",dv.cex=1.1,n=10000)
#dev.off()

# Predictions
varnames <- mod$namesx
vard <- "IndustryMachinery_Emp_PROP"
x0 <- min(X[,vard])
x1 <- max(X[,vard])
xvals <- list(c(vard,x0),c("RUS_MAJ",0))
X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
xvals <- list(c(vard,x1),c("RUS_MAJ",0))
X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
rr <- BMA.pred.alt(mod,X.0=X.0,X.1=X.1,invlink=invlink,rnd=2);rr
#print.xtable(rr,file=paste0("Output/BMA_mod1_inter_RR1a",lv,".tex"))

# Predictions
varnames <- mod$namesx
vard <- "IndustryMachinery_Emp_PROP"
x0 <- mean(X[,vard])-sd(X[,vard])
x1 <- mean(X[,vard])+sd(X[,vard])
xvals <- list(c(vard,x0),c("RUS_MAJ",0))
X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
xvals <- list(c(vard,x1),c("RUS_MAJ",0))
X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
rr <- BMA.pred.alt(mod,X.0=X.0,X.1=X.1,invlink=invlink,rnd=2);rr
#print.xtable(rr,file=paste0("Output/BMA_mod1_inter_RR1b",lv,".tex"))

# Predictions loop
varnames <- mod$namesx
vardz <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
cfmatz <- list()
for(v in 1:length(vardz)){
  vard <- vardz[v]
  x0 <- mean(X[,vard])-sd(X[,vard])
  x1 <- mean(X[,vard])+sd(X[,vard])
  cfmat <- as.data.frame(matrix(NA,4,2))
  class(cfmat[,1])
  #for(j in 1:2){cfmat[,j]<-as.character(cfmat[,j])}
  names(cfmat) <- c("Mean","CI")
  row.names(cfmat) <- c("R-X0","R-X1","X-R0","X-R1")
  cfz <- list(c(0,1,x0,x0),c(0,1,x1,x1),c(0,0,x0,x1),c(1,1,x0,x1))
  for(i in 1:length(cfz)){
    xvals <- list(c("RUS_MAJ",cfz[[i]][1]),c(vard,cfz[[i]][3]))
    X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    xvals <- list(c("RUS_MAJ",cfz[[i]][2]),c(vard,cfz[[i]][4]))
    X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    cfmat[i,1] <- as.character(BMA.pred.alt(mod=mod,X.0=X.0,X.1=X.1,invlink=invlink)[4,1])
    cfmat[i,2] <- as.character(BMA.pred.alt(mod=mod,X.0=X.0,X.1=X.1,invlink=invlink)[4,2])
  }
  cfmatz[[v]] <- cfmat
}
names(cfmatz) <- vardz
cfmatz
#print.xtable(xtable(rbind(cfmatz[[1]],rep(NA,3),cfmatz[[2]],rep(NA,3),cfmatz[[3]]),digits=2),file=paste0("Output/BMA_mod1_inter_RR2",lv,".tex"))


# Table
mm <- cbind(mod[[2]],mod[[3]][-1],mod[[4]][-1])
row.names(mm) <- mod[[1]]
row.names(mm) <- labz
mm
#print.xtable(xtable(mm,digits=2),file=paste0("Output/BMA_mod1_inter",lv,".tex"))


#####################################
# Duration until control
# (Figure 5a)
#####################################

# Merge survival times
load("Data/SurvMat_v.RData")
data <- merge(data,surv.mat,by="geonameid")
summary(data)

# BMA
X.varz <- c(ind,eth1_b,ter,dem,rus,pol,int_b)
y.varz <- c("SURV_CTR")
cens <- c("SURV_CTR_C")
datz <- na.omit(data[,c(y.varz,cens,X.varz)])
cens <- datz[,cens]
X <- datz[,X.varz]
y <- datz[row.names(datz)%in%row.names(X),y.varz]
mod <- bic.surv(x=X,surv.t=y,cens=cens,OR=10^3)
summary(mod,n.models=1,digits=1)

# Fix model weights
ix <- list(c("MCH_RUS_b","IndustryMachinery_Emp_PROP","RUS_MAJ"),c("MIN_RUS_b","IndustryMining_Emp_PROP","RUS_MAJ"),c("MTL_RUS_b","IndustryMetals_Emp_PROP","RUS_MAJ"))
mod <- bma.int(mod=mod,ix=ix); mod

# Plot
labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010","Machine x Russian","Mining x Russian","Metals x Russian")
#png(file=paste0("Output/BMA_mod3_inter",lv,".png"),width=4.5,height=4.5,units="in",res=300)
BMA.plot.surv(mod,labz,dvlab="Establishment of rebel control",dv.cex=1.1,n=10000)
#dev.off()

# Hazard ratios
varnames <- mod$namesx
vardz <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
cfmatz <- list()
for(v in 1:length(vardz)){
  vard <- vardz[v]
  x0 <- mean(X[,vard])-sd(X[,vard])
  x1 <- mean(X[,vard])+sd(X[,vard])
  cfmat <- as.data.frame(matrix(NA,4,3))
  names(cfmat) <- c("Mean","Lower","Upper")
  row.names(cfmat) <- c("R-X0","R-X1","X-R0","X-R1")
  cfz <- list(c(0,1,x0,x0),c(0,1,x1,x1),c(0,0,x0,x1),c(1,1,x0,x1))
  for(i in 1:length(cfz)){
    xvals <- list(c("RUS_MAJ",cfz[[i]][1]),c(vard,cfz[[i]][3]))
    X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    xvals <- list(c("RUS_MAJ",cfz[[i]][2]),c(vard,cfz[[i]][4]))
    X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    cfmat[i,] <- (unlist(BMA.pred.HR(mod=mod,X.0=X.0,X.1=X.1))-1)*100
  }
  cfmatz[[v]] <- cfmat
}
names(cfmatz) <- vardz
cfmatz
#print.xtable(xtable(rbind(cfmatz[[1]],rep(NA,3),cfmatz[[2]],rep(NA,3),cfmatz[[3]]),digits=2),file=paste0("Output/BMA_mod3_inter_HR",lv,".tex"))

# Hazard ratios 2
v <- 3
varnames <- mod$namesx
vardz <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
cfmatz <- list()
for(v in 1:length(vardz)){
  vard <- vardz[v]
  x0 <- mean(X[,vard])-sd(X[,vard])
  x1 <- mean(X[,vard])+sd(X[,vard])
  cfmat <- as.data.frame(matrix(NA,4,2))
  class(cfmat[,1])
  names(cfmat) <- c("Mean","CI")
  row.names(cfmat) <- c("R-X0","R-X1","X-R0","X-R1")
  cfz <- list(c(0,1,x0,x0),c(0,1,x1,x1),c(0,0,x0,x1),c(1,1,x0,x1))
  for(i in 1:length(cfz)){
    xvals <- list(c("RUS_MAJ",cfz[[i]][1]),c(vard,cfz[[i]][3]))
    X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    xvals <- list(c("RUS_MAJ",cfz[[i]][2]),c(vard,cfz[[i]][4]))
    X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    cfmat[i,1] <- unlist(BMA.pred.HR2(mod=mod,X.0=X.0,X.1=X.1))[1]
    cfmat[i,2] <- unlist(BMA.pred.HR2(mod=mod,X.0=X.0,X.1=X.1))[2]
  }
  cfmatz[[v]] <- cfmat
}
names(cfmatz) <- vardz
#print.xtable(xtable(rbind(cfmatz[[1]],rep(NA,3),cfmatz[[2]],rep(NA,3),cfmatz[[3]]),digits=2),file=paste0("Output/BMA_mod3_inter_HR2",lv,".tex"))

# Table
mm <- do.call(cbind,mod)[,-1]
mm <- as.data.frame(mm)
for(j in 1:ncol(mm)){mm[,j] <- as.numeric(as.character(mm[,j]))}
row.names(mm) <- mod[[1]]
row.names(mm) <- labz
mm
#print.xtable(xtable(mm,digits=2),file=paste0("Output/BMA_mod3_inter",lv,".tex"))


#####################################
# Loss of rebel control
# (Figure 5b)
#####################################

# Subset under rebel control
datasub <- data[data$SURV_CTR_C==1,]

# BMA
X.varz <- c(ind,eth1_b,ter,dem,rus,pol,int_b)
y.varz <- c("SURV_LIB2")
X <- na.omit(datasub[,X.varz])
y <- datasub[row.names(datasub)%in%row.names(X),y.varz]
cens <- datasub[,c("SURV_LIB2_C")]
mod <- bic.surv(x=X,surv.t=y,cens=cens,strict=F,OR=10^3)
summary(mod,n.models=1,digits=1)

# Fix model weights
ix <- list(c("MCH_RUS_b","IndustryMachinery_Emp_PROP","RUS_MAJ"),c("MIN_RUS_b","IndustryMining_Emp_PROP","RUS_MAJ"),c("MTL_RUS_b","IndustryMetals_Emp_PROP","RUS_MAJ"))
mod <- bma.int(mod=mod,ix=ix)

# Plot
labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010","Machine x Russian","Mining x Russian","Metals x Russian")
#png(file=paste0("Output/BMA_mod4_inter",lv,".png"),width=4.5,height=4.5,units="in",res=300)
BMA.plot.surv(mod,labz,dvlab="Loss of rebel control",dv.cex=1.1,n=10000)
#dev.off()

# Hazard ratios
varnames <- mod$namesx
vardz <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
cfmatz <- list()
for(v in 1:length(vardz)){
  vard <- vardz[v]
  x0 <- mean(X[,vard])-sd(X[,vard])
  x1 <- mean(X[,vard])+sd(X[,vard])
  cfmat <- as.data.frame(matrix(NA,4,3))
  names(cfmat) <- c("Mean","Lower","Upper")
  row.names(cfmat) <- c("R-X0","R-X1","X-R0","X-R1")
  cfz <- list(c(0,1,x0,x0),c(0,1,x1,x1),c(0,0,x0,x1),c(1,1,x0,x1))
  for(i in 1:length(cfz)){
    xvals <- list(c("RUS_MAJ",cfz[[i]][1]),c(vard,cfz[[i]][3]))
    X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    xvals <- list(c("RUS_MAJ",cfz[[i]][2]),c(vard,cfz[[i]][4]))
    X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    cfmat[i,] <- (unlist(BMA.pred.HR(mod=mod,X.0=X.0,X.1=X.1))-1)*100
  }
  cfmatz[[v]] <- cfmat
}
names(cfmatz) <- vardz
cfmatz
#print.xtable(xtable(rbind(cfmatz[[1]],rep(NA,3),cfmatz[[2]],rep(NA,3),cfmatz[[3]]),digits=2),file=paste0("Output/BMA_mod4_inter_HR",lv,".tex"))


# Hazard ratios 2
varnames <- mod$namesx
vardz <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
cfmatz <- list()
for(v in 1:length(vardz)){
  vard <- vardz[v]
  x0 <- mean(X[,vard])-sd(X[,vard])
  x1 <- mean(X[,vard])+sd(X[,vard])
  cfmat <- as.data.frame(matrix(NA,4,2))
  class(cfmat[,1])
  names(cfmat) <- c("Mean","CI")
  row.names(cfmat) <- c("R-X0","R-X1","X-R0","X-R1")
  cfz <- list(c(0,1,x0,x0),c(0,1,x1,x1),c(0,0,x0,x1),c(1,1,x0,x1))
  for(i in 1:length(cfz)){
    xvals <- list(c("RUS_MAJ",cfz[[i]][1]),c(vard,cfz[[i]][3]))
    X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    xvals <- list(c("RUS_MAJ",cfz[[i]][2]),c(vard,cfz[[i]][4]))
    X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    cfmat[i,1] <- unlist(BMA.pred.HR2(mod=mod,n.sim=1000,X.0=X.0,X.1=X.1))[1]
    cfmat[i,2] <- unlist(BMA.pred.HR2(mod=mod,n.sim=1000,X.0=X.0,X.1=X.1))[2]
  }
  cfmatz[[v]] <- cfmat
}
names(cfmatz) <- vardz
cfmatz
#print.xtable(xtable(rbind(cfmatz[[1]],rep(NA,3),cfmatz[[2]],rep(NA,3),cfmatz[[3]]),digits=2),file=paste0("Output/BMA_mod4_inter_HR2",lv,".tex"))

# Table
mm <- do.call(cbind,mod)[,-1]
mm <- as.data.frame(mm)
for(j in 1:ncol(mm)){mm[,j] <- as.numeric(as.character(mm[,j]))}
row.names(mm) <- mod[[1]]
row.names(mm) <- labz
mm
#print.xtable(xtable(mm,digits=2),file=paste0("Output/BMA_mod4_inter",lv,".tex"))


#####################################
# Intensity of rebel violence
# (Figure 4b)
#####################################

# Load municipality-week data
lv <- "_vw"
load(paste0("Data/Combined_1PD",lv,".RData"))
data <- data.vw
tvar <- "WID"
data <- streg.sort(data, unitvar="geonameid", timevar=tvar)
disag <- data[data[,tvar]==min(data[,tvar]),"geonameid"]

# Create lag terms
coords <- cbind(gn$longitude,gn$latitude)
IDs <- gn$geonameid
W_tri <- tri2nb(coords, row.names = IDs)
W_tri_mat <- nb2mat(W_tri, style="W", zero.policy=TRUE)
data$REB_ALL_W <- streg.slag(data, timevar=tvar, lagvar="REB_ALL", Wmat=W_tri_mat)
data$REB_ALL_t <- streg.tlag(data, timevar=tvar, lagvar="REB_ALL")
data$REB_ALL_Wt <- streg.tlag(data, timevar=tvar, lagvar="REB_ALL_W")
data$REB_ALL_bW <- streg.slag(data, timevar=tvar, lagvar="REB_ALL_b", Wmat=W_tri_mat)
data$REB_ALL_bt <- streg.tlag(data, timevar=tvar, lagvar="REB_ALL_b")
data$REB_ALL_bWt <- streg.tlag(data, timevar=tvar, lagvar="REB_ALL_bW")
data$GOV_ALL_W <- streg.slag(data, timevar=tvar, lagvar="GOV_ALL", Wmat=W_tri_mat)
data$GOV_ALL_t <- streg.tlag(data, timevar=tvar, lagvar="GOV_ALL")
data$GOV_ALL_Wt <- streg.tlag(data, timevar=tvar, lagvar="GOV_ALL_W")
data$GOV_ALL_bW <- streg.slag(data, timevar=tvar, lagvar="GOV_ALL_b", Wmat=W_tri_mat)
data$GOV_ALL_bt <- streg.tlag(data, timevar=tvar, lagvar="GOV_ALL_b")
data$GOV_ALL_bWt <- streg.tlag(data, timevar=tvar, lagvar="GOV_ALL_bW")
data$CONTROL_ANY <- 1*(data$CONTROL_MEAN>0)
data$CONTROL_W <- streg.slag(data, timevar=tvar, lagvar="CONTROL_ANY", Wmat=W_tri_mat)
data$CONTROL_t <- streg.tlag(data, timevar=tvar, lagvar="CONTROL_ANY")
data$CONTROL_Wt <- streg.tlag(data, timevar=tvar, lagvar="CONTROL_W")
rm(data.vw,W_tri,W_tri_mat)

# Rescale variables
data$POP <- data$POP/max(data$POP,na.rm=T)
data$dem <- data$dem/max(data$dem,na.rm=T)
data$DIST2ROAD <- data$DIST2ROAD/max(data$DIST2ROAD,na.rm=T)
data$DIST2RUS <- data$DIST2RUS/max(data$DIST2RUS,na.rm=T)
data$RUS_CENSUS2001 <- data$RUS_CENSUS2001/100
data$V2010A_YANUKOVYCH <- (data[,"V2010A_YANUKOVYCH"]-min(data[,"V2010A_YANUKOVYCH"]))/(max(data[,"V2010A_YANUKOVYCH"])-min(data[,"V2010A_YANUKOVYCH"]))

# Create interactions
data$RUS_MAJ <- 1*(data$RUS_CENSUS2001>.5)
data$MCH_RUS <- data$IndustryMachinery_Emp_PROP*data$RUS_CENSUS2001
data$MIN_RUS <- data$IndustryMining_Emp_PROP*data$RUS_CENSUS2001
data$MTL_RUS <- data$IndustryMetals_Emp_PROP*data$RUS_CENSUS2001 
data$MCH_RUS_b <- data$IndustryMachinery_Emp_PROP*data$RUS_MAJ
data$MIN_RUS_b <- data$IndustryMining_Emp_PROP*data$RUS_MAJ
data$MTL_RUS_b <- data$IndustryMetals_Emp_PROP*data$RUS_MAJ 

# Choose covariates
lags <- c("REB_ALL_Wt","REB_ALL_t")
eth1 <- c("RUS_CENSUS2001")
eth1_b <- c("RUS_MAJ")
int <- c("MCH_RUS","MIN_RUS","MTL_RUS")
int_b <- c("MCH_RUS_b","MIN_RUS_b","MTL_RUS_b")

# BMA
X.varz <- c(ind,eth1_b,ter,dem,rus,pol,int_b,lags)
y.varz <- c("REB_ALL")
X <- na.omit(data[,X.varz])
y <- data[row.names(data)%in%row.names(X),y.varz]
mod <- bic.glm(x=X,y=y,glm.family="quasipoisson",strict=F,OR=10^4)
summary(mod,n.models=1,digits=1)
invlink <- mod$linkinv

# Fix model weights
ix <- list(c("MCH_RUS_b","IndustryMachinery_Emp_PROP","RUS_MAJ"),c("MIN_RUS_b","IndustryMining_Emp_PROP","RUS_MAJ"),c("MTL_RUS_b","IndustryMetals_Emp_PROP","RUS_MAJ"))
mod <- bma.int(mod=mod,ix=ix)

# Predictions
vard <- "IndustryMachinery_Emp_PROP"
x0 <- 0
x1 <- 1
xvals <- list(c(vard,x0),c("RUS_MAJ",1))
X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
xvals <- list(c(vard,x1),c("RUS_MAJ",1))
X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
rr <- BMA.pred.alt(mod,X.0=X.0,X.1=X.1,invlink=invlink,rnd=2);rr
#print.xtable(rr,file=paste0("Output/BMA_mod2_inter_RR1a",lv,".tex"))

# Predictions
vard <- "IndustryMachinery_Emp_PROP"
x0 <- min(X[,vard])
x1 <- max(X[,vard])
xvals <- list(c(vard,x0),c("RUS_MAJ",1))
X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
xvals <- list(c(vard,x1),c("RUS_MAJ",1))
X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
rr <- BMA.pred.alt(mod,X.0=X.0,X.1=X.1,invlink=invlink,rnd=2);rr
rr
#print.xtable(rr,file=paste0("Output/BMA_mod2_inter_RR1b",lv,".tex"))

# Predictions
vard <- "IndustryMachinery_Emp_PROP"
x0 <- mean(X[,vard])-sd(X[,vard])
x1 <- mean(X[,vard])+sd(X[,vard])
xvals <- list(c(vard,x0),c("RUS_MAJ",1))
X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
xvals <- list(c(vard,x1),c("RUS_MAJ",1))
X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
rr <- BMA.pred.alt(mod,X.0=X.0,X.1=X.1,invlink=invlink,rnd=2);rr
#print.xtable(rr,file=paste0("Output/BMA_mod2_inter_RR1c",lv,".tex"))

# Predictions loop
varnames <- mod$namesx
vardz <- c("IndustryMachinery_Emp_PROP","IndustryMining_Emp_PROP","IndustryMetals_Emp_PROP")
cfmatz <- list()
for(v in 1:length(vardz)){
  vard <- vardz[v]
  x0 <- mean(X[,vard])-sd(X[,vard])
  x1 <- mean(X[,vard])+sd(X[,vard])
  cfmat <- as.data.frame(matrix(NA,4,2))
  class(cfmat[,1])
  names(cfmat) <- c("Mean","CI")
  row.names(cfmat) <- c("R-X0","R-X1","X-R0","X-R1")
  cfz <- list(c(0,1,x0,x0),c(0,1,x1,x1),c(0,0,x0,x1),c(1,1,x0,x1))
  for(i in 1:length(cfz)){
    xvals <- list(c("RUS_MAJ",cfz[[i]][1]),c(vard,cfz[[i]][3]))
    X.0 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    xvals <- list(c("RUS_MAJ",cfz[[i]][2]),c(vard,cfz[[i]][4]))
    X.1 <- BMA.X(X=X,xvals=xvals,varnames=varnames,ix=ix)
    cfmat[i,1] <- as.character(BMA.pred.alt(mod=mod,X.0=X.0,X.1=X.1,invlink=invlink)[4,1])
    cfmat[i,2] <- as.character(BMA.pred.alt(mod=mod,X.0=X.0,X.1=X.1,invlink=invlink)[4,2])
  }
  cfmatz[[v]] <- cfmat
}
names(cfmatz) <- vardz
cfmatz
#print.xtable(xtable(rbind(cfmatz[[1]],rep(NA,3),cfmatz[[2]],rep(NA,3),cfmatz[[3]]),digits=2),file=paste0("Output/BMA_mod2_inter_RR2",lv,".tex"))


# Plot
labz <- c("Machine industry","Mining industry","Metals industry","Russian language","Forest","Elevation","Distance to road","Population density","Distance to Russia","Yanukovych vote, 2010","Machine x Russian","Mining x Russian","Metals x Russian",expression(paste(bold(W)," Rebel Violence",s[t-1])),expression(paste("Rebel Violence",s[t-1])))
#png(file=paste0("Output/BMA_mod2_inter",lv,".png"),width=4.5,height=4.5,units="in",res=300)
BMA.plot(mod,labz,dvlab="Intensity of rebel violence (weekly)",dv.cex=1.1,n=10000)
#dev.off()

# Table
mm <- cbind(mod[[2]],mod[[3]][-1],mod[[4]][-1])
row.names(mm) <- mod[[1]]
row.names(mm) <- labz
#print.xtable(xtable(mm,digits=2),file=paste0("Output/BMA_mod2_inter",lv,".tex"))
